Evolving black hole-neutron star binaries in general relativity using pseudospectral 

and finite difference methods 
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We present a code for solving the coupled Einstein-hydrodynamics equations to evolve relativistic, 
self-gravitating fluids. The Einstein field equations are solved in generalized harmonic coordinates 
on one grid using pseudospectral methods, while the fluids are evolved on another grid using shock- 
capturing finite difference or finite volume techniques. We show that the code accurately evolves 
equilibrium stars and accretion flows. Then we simulate an equal-mass nonspinning black hole- 
neutron star binary, evolving through the final four orbits of inspiral, through the merger, to the 
final stationary black hole. The gravitational waveform can be reliably extracted from the simulation. 
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I. INTRODUCTION 

Compact object binaries containing neutron stars [i.e. 
neutron star-neutron star (NSNS) and black hole-neutron 
star (BHNS) binaries] are perhaps as important in mod- 
ern astrophysics as binary black holes. Both BHNS and 
NSNS binaries should be excellent sources of gravita- 
tional waves for the ground-based interferometric detec- 
tors LIGO, GEO, VIRGO, and TAMA. It is possible [J 
that the detection rate of these binaries will actually 
be greater than that of binary black holes. NSNS and 
BHNS binaries are also interesting because they are lead- 
ing candidates for explaining the production of short- 
duration gamma-ray bursts (GRBs) [2j, especially given 
observations that locate some short GRBs in elliptical 
galaxies 0, 0] or rule out an associated supernova ||. 
NSNS and BHNS mergers may also be important for un- 
derstanding the observed abundances of the heavy ele- 
ments that are formed by rapid neutron capture in the 
r- process @. 

NSNS and BHNS binary mergers can be accurately 
modeled only by numerical simulations. In such systems, 
the spacetime metric and neutron star fluid are both dy- 
namical. They are also strongly coupled, and hence must 
be evolved simultaneously. In addition, it is clear from 
the high compactions of the binary objects that only sim- 
ulations in full general relativity will be adequate. For 
NSNS binaries, such simulations are particularly needed 
to study the possible collapse of the post-merger rem- 
nant. For BHNS binaries, the crucial questions are 1) 
whether or not the neutron star is tidally disrupted be- 
fore it plunges into the black hole, and 2) if the neutron 
star does disrupt, what fraction of the star's matter is 
promptly swallowed by the hole, what fraction is ejected, 
and what fraction forms an accretion disk. To produce a 
GRB, a substantial accretion disk must remain. To con- 
tribute to the abundance of r-process elements, matter 
must be expelled from the system. Studies strongly sug- 
gest that the non-Newtonian form of the gravitational po- 
tential 0] , gravitational radiation reaction Q , and black 



hole spin Q all significantly affect the size of the post- 
merger accretion disk, underscoring the need for fully 
relativistic simulations. 

Of the two types of systems, NSNS binaries are better 
studied. Newtonian simulations have been performed us- 
ing realistic equations of state together with neutrino ra- 
diation effects [13, EH, EEHIl and magnetic fields [H[. 
Simulations using complicated equations of state have 
also been carried out using the conformally flat approx- 
imation to general relativity [l6l |. General relativistic 
simulations have been performed using T-law equations 
of state [HI GJ, [IlJSSJp, HI and using more realistic 
equations of state 23, 24|. Recently, general relativistic 
merger evolutions have been performed that include the 
neutron star magnetic field [SB], [26J . 

Numerical modeling of BHNS binary mergers has been 
carried out using Newtonian or pseudo-Newtonian grav- 
ity (in which the black hole is represented by a point 
mass) using T-law [13, HI] and realistic nuclear 0, [H, [3(| 
equations of state. BHNS simulations have also been 
done in conformal gravity [3l|, although so far only in 
the extreme mass ratio limit, in which the black hole 
is much more massive than the neutron star. The first 
fully relativistic BHNS simulation was of a head-on colli- 
sion [32j]. Most recently, two groups have independently 
evolved configurations of BHNS binaries in full general 
relativity from quasi-circular inspiral through merger. 
One group includes Shibata, Taniguchi, Uryii, and Ya- 
mamoto [H, [13, The members of the other group 
are Etienne, Faber, Liu, Shapiro, Taniguchi, and Baum- 
garte (36[ . The evolutions produced by these two groups 
agree qualitatively, with both now finding very small 
post-merger disks. However, both groups have so far 
restricted themselves to V = 2 polytropic equations of 
state and to initially non-spinning black holes. The case 
of spinning black holes has been studied by Rantsiou, 
Kobayashi, and Laguna Q using a Kerr background met- 
ric and a Newtonian approximation for the neutron star 
self-gravity. They find that black hole spin has a strong 
influence on the size of the post-merger accretion disk. 
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For both types of binary, but especially BHNS binaries, 
the parameter space remains poorly explored. Also, there 
is an unmet need for results to be checked by multiple 
independent codes. 

In each of the above-mentioned calculations in which 
the metric variables are truly dynamical fields, these 
fields were evolved numerically using finite differencing 
(FD). For a stable evolution, such algorithms should con- 
verge to the exact solution as some power of the grid 
spacing. Also, shock-capturing techniques have been de- 
veloped which allow FD codes to evolve fluids with dis- 
continuities stably and accurately. FD codes usually re- 
quire very large grids in order to obtain accurate results, 
although this problem can be mitigated by using mesh 
refinement or high-order differencing. 

Einstein's equations can also be evolved using spectral 
methods. In this case, functions are approximated as 
truncated series expansions in a set of orthogonal basis 
functions. Derivatives of the approximated functions are 
then computed exactly. The pseudospectral (PS) method 
is a type of spectral method that uses the values of func- 
tions on a spatial grid of collocation points, rather than 
directly using the spectral coefficients. This has the ad- 
vantage that pointwise operations are as straightforward 
in PS methods as they are in FD methods. In multido- 
main PS methods, the computational region is divided 
into domains, each with its own basis functions and corre- 
sponding collocation points. For smooth functions, spec- 
tral methods (including PS) converge exponentially to the 
exact solution as the number of basis functions (or, for PS 
methods, the number of collocation points) is increased. 
This allows PS methods to get accurate results with much 
smaller grids than those used by FD codes. A PS code 
for solving the Einstein equations has been developed by 
the Cornell-Caltech relativity group [H H, IM S3- lt 
has been used to simulate the inspiral of binary black 
holes for many orbits with v ery high accuracy for a fairly 
low computational cost 0, SH 53 ■ 

There is a difficulty, however, in extending PS meth- 
ods to evolving non-vacuum spactimes. Because of the 
presence of stellar surfaces and, in some cases, hydro- 
dynamic shocks, the evolved variables are not smooth 
in all derivatives. In these cases, spectral representa- 
tions do not converge exponentially to the exact solution. 
Rather, they display Gibbs oscillations near the discon- 
tinuity that converge away only like some power of the 
number of collocation points, with the order of conver- 
gence given by the order of the discontinuity. The oscil- 
lations can be controlled using special forms of filtering 
or artificial viscosity (e.g. SaSlI) but exponential con- 
vergence is still lost. In some cases, the problem can be 
avoided by placing domain boundaries at discontinuities, 
e.g. at the stellar surface. However, this is not prac- 
tical for complicated shocks or when stars become very 
deformed (e.g. during tidal disruption). 

In this paper, we use the "mixed" approach developed 
by Dimmelmeier et al [46j for a conformal gravity code 
which has been used quite successfully to study super- 



nova core collapse [13, IH, [49|, [5(| , and extend it to full 
general relativity. In this method, the spacetime metric 
is evolved on one grid using PS methods, the fluid vari- 
ables are evolved on a second grid using shock-capturing 
FD or finite volume techniques, and the two grids com- 
municate by interpolation. Below, these two grids are 
referred to as the "PS grid" and the "fluid grid" . 

This approach would seem to utilize the strongest fea- 
tures of each method. The spacetime is expected to be 
smoother than the fluid variables (e.g. at a stellar sur- 
face, the discontinuity in metric components appears at 
a higher derivative than in the density), and so PS tech- 
niques should work better on the metric than the fluid. 
Also, our code uses many spectral domains, and we ex- 
pect discontinuities to appear in only a few domains. In 
the domains without discontinuities, the functions can 
be represented spectrally with the same accuracy as in 
smooth problems. In the domains with discontinuities, 
convergence will be limited to a power law. Of course, 
this error will propagate to the other domains. However, 
the domain decomposition can often be chosen so that 
the slower converging domains take up a small fraction 
of the overall computational region, and their effect on 
the overall accuracy is correspondingly small. We can 
also use higher resolution in these domains if the error 
is still too large. Unlike the strategy of fitting domain 
boundaries to discontinuities, we do not require the ex- 
act locations of discontinuities, but only approximate lo- 
cations. 

Even if the evolution of the fields converges rapidly, 
spectral accuracy is still lost because of the evolution of 
the fluids, which will at best converge as a fixed power in 
the grid spacing on the fluid grid. However, this is not 
a problem if the resolution on this grid is high. And, in 
fact, our mixed technique algorithm allows us to achieve 
high resolution in the hydrodynamic evolution at a sur- 
prisingly low computational cost. This is because the 
fluid grid only needs to cover the region containing the 
matter. This can provide a huge savings for binary inspi- 
ral calculations. For a BHNS inspiral, the fluid grid can 
be a box centered on the neutron star with outer bound- 
aries slightly outside the star. For a NSNS inspiral, one 
would need two such boxes. For a BHNS system in which 
the star is shedding mass onto the black hole, the fluid 
grid would have to cover a region containing both binary 
objects, but even then, it would not have to extend all 
the way out into the wave zone, as it would if the metric 
were being evolved on the same grid. 

In this paper, we describe our code and test its abil- 
ity to evolve BHNS binaries. In Section HH we describe 
our evolution code. In Section IIII1 we present tests of 
this code. Next, we apply our code to model the inspi- 
ral and merger of a BHNS binary. We evolve an equal 
mass binary with an initially non-spinning black hole and 
irrotational neutron star. Section HVl describes the inspi- 
ral calculation, with particular emphasis on the accuracy 
and rate of convergence of these simulations. Section [V] 
describes the merger. Finally, Section IVII gives our con- 
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elusions and future directions for our work. 



II. EVOLUTION CODE 



A. Evolution of the Spacetime 



We evolve Einstein's 
harmonic formulation 




uations using the generalized 
We use a first order rep- 
resentation of the system [39( , in which the fundamental 
variables are the spacetime metric ip a bi its spatial first 
derivatives <5>i a b, and its first derivatives in the direction 
normal to the slice Yl a b- (Throughout this paper, Latin 
indices from the first part of the alphabet a, b, c, . . . run 
from to 3, while k, . . . run from 1 to 3.) From these 
functions, one can easily extract the 3- metric gij, shift 
f3 l , lapse a, and extrinsic curvature ify. The evolution 
of the gauge is determined by the gauge source functions 
H a = —ip cd r aci i, which are freely specifiable functions of 
space and time. The evolution equations for ip a b, <&i a b, 
and n o b are as given in our earlier paper [39[, except we 
add the matter source term 



d t U ab = 2a(T ab - -ipabT^a) 



(1) 



We evolve ip ao , ^iab, and n a b using a multi-domain 
pseudospectral (PS) code described in earlier papers [13, 
|40L l4l| . Each PS domain is either a spherical shell, a 
cylindrical shell, a cube, a "cubed sphere" [53[, a filled 
sphere, or a filled cylinder. Spherical harmonics Yi m are 
used as angular basis functions on spheres. Fourier func- 
tions e lm ^ are used for the azimuthal direction on cylin- 
ders. Chebyshev polynomials T n are used as basis func- 
tions for each direction on cubes and cubed spheres, for 
the radial direction on spherical shells, and for the radial 
and z-directions on cylindrical shells. The basis functions 
on filled spheres and cylinders must be chosen specially 
to have the proper behavior at the origin and axis, re- 
spectively. For example, for spheres, we must decompose 
a function / as f(r,d,(/>) = J2 nlrn Q n i{r)Y hn (6, (f) (note 
the coupled indices), where Q n i cx r l near the origin. We 
use the functions introduced by Matsushima and Mar- 
cus [HI with a = 1, /? = 1 for filled cylinders and a = 1, 
(3 = 2 for filled spheres. We use boundary conditions [HH 
designed to prevent the influx of gravitational radiation 
and constraint violating metric perturbations. Our do- 
main decomposition is chosen to leave an unfilled (ex- 
cised) region inside the black hole. No explicit boundary 
condition need be applied at the boundary of the excised 
region because all characteristics there flow out of the 
grid. 

The PS grid supplies the fluid grid with the metric 
fields a, (3 Z , g^, Kij, dice, Sift, and dig^ k . These func- 
tions must be interpolated from the PS gridpoints to the 
fluid gridpoints. This could be done directly at each fluid 
gridpoint by summing the values of all the basis functions 
at that point with weights given by the known spectral 



coefficients. This would, however, be prohibitively ex- 
pensive, as it would involve ~ NpjyNpg operations, with 
Nfd the number of destination points (i.e. the number 
of points on the fluid grid) and Nps the number of points 
on a domain of the PS grid. Instead, we used a trick in- 
troduced by Boyd (56|. We first interpolate the metric 
fields onto a finer PS grid. (We usually triple the num- 
ber of collocation points.) This can be done cheaply by 
switching to spectral space (often an Nps log Nps oper- 
ation procedure), adding more basis functions with zero 
coefficients, and switching back to physical space. Then, 
one can get an accurate estimate of the field at fluid grid- 
points by doing polynomial interpolation from the refined 
PS grid, which takes only ~ Npp> operations. Using this 
procedure, the CPU time spent on interpolation, while 
still significant, is smaller than the time spent on evolv- 
ing on the PS and fluid grids. 



B. Evolution of the hydrodynamic fields 

We model the neutron star matter as a perfect fluid 
with rest mass density p, pressure P, specific internal 
energy e, and 4- velocity u a , so that the stress tensor is 



T ab = phu a U b + Pijj ab 



(2) 



where h = 1 + e + P/p is the specific enthalpy. The evo- 
lution of the fluid is determined by the laws of baryon 
conservation (pu a )- a = and energy-momentum conser- 
vation T ab [a — 0. These give five evolution equations 
for the variables D — ay/gu°p, t — a 2 ^/gT 00 — D, and 
Sk = a^/gT°k- Here u° is given by the normalization 
condition ipabU a u b = — 1. The pressure is given by the 
equation of state. Our code is written to evolve general 
equations of state of the form 



e = £coid (p) + eth 
P = Pcoid(p) + (r th - l)pe th 



(3) 
(4) 



(c.f. [23j|). In this paper, we will mainly use the simple 
T-law equation of state: 



Pcoid = np r 

£cow = P C oid/[p(r - 1)] 
r th = r , 



(5) 
(G) 
(7) 



which is equivalent to 

P = (r - l)pe . (8) 

The hydrodynamic equations have the form 

d t u + d z F l =S, (9) 

where u is called a "conservative variable" (D, r, or Sk, 
in our case), while F l and S are the associated flux and 
source terms, respectively. To solve these equations, one 
first divides the computational domain into cells, with 
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one cell associated with each gridpoint. Then, from the 
values of u in the cells, one determines the fluxes F l 
at the interfaces between cells. From these, one can 
compute the net flux into each cell, hence obtaining 
the diF 1 term in Eq. ([9]). A specific hydrodynamics 
evolution algorithm is determined by a set of choices: 

1) Discretization 

The values of u at gridpoints can be chosen to be 
the values of u at the centers of cells {finite difference 
discretizaton) or to be the average values of u inside the 
cells {finite volume discretizaton). The distinction only 
matters for schemes that are higher than second-order 
accurate. All of our higher order schemes assume finite 
differencing except for PPM (see below), for which we 
have coded finite difference and finite volume options. 
Another discretization issue is the shape of the grid. We 
always use uniformly spaced Cartesian meshes, i.e., at 
grid indicies (i, j,k) we have gridpoints (cell centers) 
Xi t j tk = (iAx + x ,jAy + y , kAz + z ). 

2) Interpolation variables 

To compute the fluxes, five independent variables must 
be interpolated from cell centers or averages to cell faces. 
We interpolate p to get face densities and Ui to get in- 
formation on face velocities. We have experimented with 
several choices for the fifth variable, which should carry 
information about face internal energies. The pressure 
P is a common choice, and we have found that it works 
adequately. However, it may not be ideal for evolving 
low-temperature stars, the main application of this pa- 
per. This is because interpolating P does not separate 
the zero-temperature component of P from the thermal 
component. So, even if P = P co \d{p) exactly on all grid 
values, the interpolated P and p will not satisfy this re- 
lationship because of interpolation error, i.e., there is in- 
terpolation "heating" . Hence, neutron star simulations 
carried out using an isentropic equation of state {P = np T 
at all times) can be very accurate for certain problems 
(e.g. [57], The isentropic treatment completely re- 

moves spurious heating, but it is not valid in the presence 
of shocks. 

Another choice would be to evolve k = P/p r - Since 
this is a constant for cold stars, it can be interpolated 
exactly. The pressure at faces is then constructed from 
the interpolated k and interpolated p. We find that 
interpolating k works well for cold stars, but behaves 
badly at shocks. This can be fixed by switching the 
scheme at discontinuities back to interpolating P. Even 
so, its usefulness is limited to polytropes. An even better 
choice is to interpolate Pth = P— P co \d, the thermal part 
of the pressure. Then the pressure on the face is the sum 
of the interpolated P t h and the P C oid computed from 
the interpolated p. With this choice, interpolation error 
introduces no heat on the faces for a zero-temperature 
star. In addition, it seems to have no problems evolving 
shocks, and it can be used with any equation of state. 
By interpolating either k or P t h, we find that spurious 



heating in our stars is often significantly reduced. 

3) Interpolation method 

Interpolation must be done in a special way to accom- 
modate the possibility of shocks. (The process is usually 
called "reconstruction" in the literature.) Consider a 
one-dimensional problem, in which we know grid values 
Pi of the function p{x). The reconstruction step involves 
computing p L = p;_i/ 2 _ e and p R = Pi-i/ 2+i , i.e. the 
values of p to the left and right of the grid cell interface. 
(In three dimensions, one must reconstruct at faces in 
each direction.) Many techniques have been developed 
for interpolating from cell centers or averages to faces 
We have implemented second-order monotonized central 
(MC) reconstruction, third-order piecewise parabolic 
(PPM) [59( reconstruction, third-order convex essentially 
nonoscillatory (CENO) [60| reconstruction, and third- 
order weighted essentially nonoscillatory (WENO) [6l[ 
reconstruction. We have found that finite volume PPM 
usually gives the best results when evolving equilibrium 
stars. We have checked that this scheme achieves third- 
order convergence for one-dimensional smooth flows. 
For more general problems, we usually find second-order 
convergence. (See tests below.) 

4) Flux computation 

From pr and pl, one can compute ur and ul, the face 
values of the conservative variables, and Fr and Fl, the 
face values of the fluxes. The "true" flux at the interface 
must be constructed by some sort of combination of the 
values on the two sides. We use the scheme of Harten, 
Lax, and van Leer (HLL) [f^l, in which the flux is 

t-! ^min-^i?, ^-max-^L CminCmax(^i?. ^l) / in \ 

*<-l/2 = . , (10) 

wnin t ^ max 

where c m ; n and c max are the left-going and right-going 
sound speeds. The first two terms in Eq. (|f Op provide an 
average value of the flux, while the last two contribute to 
a stabilizing diffusive term. 

The time derivative of the variable u is equal to 
S — divF (c.f. Eq. ©). For a finite volume method, 
the flux divergence is just 

Aa; 

plus the equivalent terms for the other two dimensions. 
For a finite difference method, other terms must be added 
to achieve third-order accuracy (see e.g. [H| for details). 

The vacuum region outside the star or stars requires 
special treatment. We fill the vacuum with a very low 
density ( 10~ 7 of the maximum value of p) atmosphere, as 
is usually done for these types of problems. At each step, 
we apply a floor on p. We also apply both a floor and a 
ceiling on the temperature in the atmosphere, requiring 
< T < 8(P co id/p). Finally, we limit the velocity in the 
atmosphere to some fraction (e.g. 0.8) of the speed of 
light. The pressure and velocity cutoffs are only applied 
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in regions where the density is below a chosen threshold 
(usually ~ 1CP 3 of the maximum value of p). Therefore, 
they should have little effect on evolutions. 

The PS grid requires T a b at each collocation point, so 
we must interpolate from the fluid grid to the PS grid. 
We have written an arbitrary-order polynomial interpo- 
lator and a third-order WENO interpolator. One would 
expect the latter to give more sensible results near strong 
shocks, but for the applications reported here, both work 
sufficiently well. 

When evolving black hole systems, we excise the grid- 
points in a neighborhood of the singularity and evolve the 
region near the excision zone using one-sided differencing 
as described in 64]. Of course, the excised region of the 
fluid grid must match that of the PS grid. Also, when 
interpolating from the fluid to the PS grid, we choose 
interpolation stencils which avoid the excised region. 

When symmetries are present in a problem, they can 
sometimes be used to reduce the size of the fluid grid. 
We have coded options to impose reflection symmetries 
on the x = 0, y = 0, or z = planes, or any combination 
of them. In these cases, the fluid grid covers only part of 
the system, with fluid quantities in other regions being 
filled for the PS grid using the appropriate symmetries. 
We have also coded a version of the hydrodynamics code 
that assumes axisymmetry and uses a two-dimensional 
fluid grid. 

Both the PS and the hydrodynamics modules of our 
code use multiple domains. This allows the code to be 
run in parallel by assigning different domains to different 
processors. We place ghost zones at processor interfaces 
in the fluid grid, so that numerical derivatives can be 
correctly taken at all true grid points. 



where (x c , y c , 0) is the point at which x l = x l , which we 
choose to be fixed for all time, while $ and a are functions 
of time. One evolves $ to correct for the binary's orbital 
rotation, and a to correct for any change in the orbital 
separation (i.e. inspiral or eccentricity). This is done by 
introducing a feedback mechanism into the evolution of 
$ and a (see [4l[ for details). 

The dual frame mappings actually used in our simula- 
tions are very similar to Eq. (|12p with one notable differ- 
ence. For large outer radius Rq, small changes in a lead 
to large absolute changes in ai?o, causing an instability. 
Therefore, we replace the linear scaling 7 = ar used in 
Eq. (TT2")) with T = ar + (1 — a)r 3 / Rq, so that we recover 
the linear scaling when r <C Rq, but the outer boundary 
is fixed even when a changes. Such a mapping has been 
successfully used for long-inspiral binary black hole sim- 
ulations [431 ] . Hereafter, we will refer to this coordinate 
map, a rotation combined with a nonlinear scaling, as 
US. 

Since the equations we integrate are written in the in- 
ertial frame, they involve derivatives with respect to x l . 
These derivatives are obtained by first taking derivatives 
with respect to grid coordinates x % and then transform- 
ing them using a Jacobian matrix. So, for example the 
equation 



dtu a + A ka pdj:u b = 



becomes 
d t vF 



dx l . 9x^ a1 ^ 



dt 



-S a -, 



dx k 



div? = 



(13) 



(14) 



Both the metric and the fluid can be evolved in this way. 
Thus Eq. ((9]) becomes 



C. Evolution of the coordinates 

Accuracy in a multidomain PS code strongly depends 
on having a domain decomposition which matches the 
"shape" of the fields being evolved. Also, the excised 
region of a black hole must remain inside the horizon. 
Thus, it is desirable for the grid to move with the black 
holes and neutron stars. We do this using the "dual co- 
ordinate frames" system developed for binary black hole 
evolutions [41| . The coordinate frame x l is set to be 
an asymptotically flat, inertial frame. All tensor compo- 
nents are evaluated with respect to this frame. The PS 
and fluid gridpoints are fixed in the computational frame 
x l . By means of a mapping between the frames, the com- 
putational coordinates approximately co-move with the 
system. For example, one can track a binary using a 
simple combination of rotation and radial scaling: 



t = t 

x = x c + a[(x - x c ) cos($) - (y - y c ) sin($)] 
V = y c + a[(x - x c ) sin($) + (y - y c ) cos($)] 
Z = az , 



(12) 



3x l 
dx k 



u di 



dx l dx k 
dt dx l 

' dx l dx k 
dt dx* 



F 



(15) 



S 



Note that we have moved the coordinate advection term 
^fr^^u into the derivative containing the flux so it can 

dx' dt ° _ 

mostly cancel the velocity advection term in F k . 

There are some subtleties involved in applying dual 
frames to our shock-capturing techniques; a straightfor- 
ward evolution of Eq. (fT5|) will often be unstable. How- 
ever, we have implemented two techniques for evolving 
fluids in moving coordinates that have proven to be ro- 
bustly stable. We believe that the advantages of moving 
coordinates are so great that it is useful to report both 
methods. 

The conceptually simpler technique is to transform the 
metric and the field variables into moving coordinates, 
and then evolve as usual in this coordinate system. One 
must remember that D, r, and Su are densities, so trans- 
forming them to the moving frame involves multiplying 



by det(dx l /dx k ). Since the spectral code evolves tensors 
using inertial frame components, coordinate transforma- 
tions must be performed when the two grids communi- 
cate. 

It is also possible to evolve the hydrodynamic equa- 
tions with inertial frame components using Eq. (|15p . 
However, one must make a few adjustments to the shock- 
capturing code. 1) The c min and c max in Eq. (jTUJ) should 
be computed in the moving frame. 2) The fluxes in 
Eq. (|10[) must be broken up into two pieces. For example, 
for the x-interface fluxes 



1.5 



•t+l/2 
1+1/2 



c) 1 (Cmin^ 



mm^m ax 



+ c 



UL) 



(16) 
(17) 



The k index may seem unfamiliar; it is not usually needed 
because in evolutions without dual frames, one only needs 
the x-derivative of F x , the y-derivative of F y , and the z- 
derivative of F z . However, the Jacobian in Eq. (fTS"]) mixes 
derivatives, so the other derivatives are now needed. To 
the two flux pieces correspond two pieces of the flux 
derivatives 



A*" 1 (^1)2 



pfc(l) V 

r i-l/2> 



(18) 
(19) 



and similarly for y and z. Next, the Jacobian in Eq. (|15p 
is applied to (d t F k )^ to give (djF k )^ . The Jacobian 
is not applied to (diF) 1 ^; this term is added directly to 
the time derivative of u. If the Jacobian is applied to 
(diFy( 2 \ it will not behave correctly as a diffusion term, 
and the code will be unstable. 



III. CODE TESTS 
A. Shocks 

As a first test of the hydrodynamics code, we perform 
a standard Riemann shock tube problem. We evolve a 
r = 5/3 fluid with an interface separating regions of dif- 
ferent p and P. The initial state is (p, P, v) — (10, 13.3, 0) 
on one side of the interface and (p,P,v) — (1,0,0) on the 
other. The interface is chosen to lie along the xy diago- 
nal, so the problem is numerically two-dimensional. We 
evolve using the PPM reconstructor on a grid of 150 2 
points and grid spacing Aa; = 0.013. The results at 
t = 0.35 are shown in Fig. [T] We see that our code 
captures all of the features of the exact solution. 



B. Spherical Accretion 

Next, we check that our code can accurately model 
accretion onto a nonrotating black hole. Note that ac- 
cretion test problems provide very sensitive checks of a 
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FIG. 1: A snapshot at t = 0.35 of the density, pressure, and 
3-velocity at each gridpoint on a line perpendicular to the 
shock interface. The lines show the exact solution. Circles, 
squares and triangles show the numerical pressure, density, 
and velocity, respectively. 



code's treatment of the effects of curved spacetime on 
the fluid flow (the non-Minkowski pieces of the fluid equa- 
tions), and also checks the code's ability to advect matter 
through an excision boundary without losing stability or 
accuracy. 

We perform three types of accretion tests. For the 
first test, we check that the code can maintain an equi- 
librium radial accretion flow. As an initial state, we use 
the well-known exact solution of Michel 65]. We use 
Kerr-Schild coordinates and excise all gridpoints inside a 
radius of r = 1.6M, where M is the mass of the black 
hole, so the event horizon is at r = 2M. We choose a 
r = 5/3 fluid and a flow with the sonic point at r s = 8M. 
As an outer boundary condition, wc hold the fluid vari- 
ables at r > 12M to the exact values. The metric fields 
are held fixed — otherwise, the black hole would grow, 
and the problem would not be stationary. Derivatives 
of the metric are computed analytically. The fluid is 
evolved to t = 100M on three grids, with grid spacings 
of Ax = 0.3M, Aa; = 0.2M, and Aa; = 0.15M. We as- 
sume that octant symmetry is maintained, and so evolve 
only one octant. The results are shown in the top panel 
of Fig O The solution remains accurate everywhere, in- 
cluding near the sonic point. We find slightly worse than 
second-order convergence. This is because of the non- 
smooth flow at the sonic point. We find that in most 
of the smooth regions, the flow converges to the exact 
solution at second order. 

For the second test problem, we perturb the initial 
state by multiplying the density by \-e~( T /( 6M >' . We 
then evolve with the same outer boundary condition as 
before. The flow should relax to the known equilibrium 
solution, and one can see from the middle panel of Fig [2] 
that it does. 

Finally, we simulate a case in which the black hole is 
moving relative to the fluid at large distance. We choose 
the fluid to have the equation of state P = p. For this 
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FIG. 2: The deviation of the density p from its exact 
equilibrium p cx . Ldcv(/) is a relative L2 norm: Ldev(/) = 
L2(/ — / cx )/L2(/ ex ), where the norms sum only over points 
outside the horizon. The top panel is a relativistic Bondi flow, 
the middle panel is a perturbed Bondi flow, and the bottom 
panel is a moving black hole flow. The deviations in the top 
and bottom panels are scaled for second-order convergence. 



equation of state, the equilibrium has been computed 
analytically by Petrich, Shapiro, and Teukolsky [66| . We 
take the relative velocity between the hole and the dis- 
tant fluid to be 0.6c in the z— direction, and we evolve in 
the black hole's frame. The fluid is evolved at two reso- 
lutions, corresponding to Ax = 0.3M and Ax = 0.2M. 
Because of the front-back asymmetry, we must evolve 
two octants. The result is shown in the bottom panel of 
Figdl The convergence is second-order. We note that for 
this equation of state, the speed of sound is equal to the 
speed of light, so there is no sonic point, and the solution 
is smooth everywhere. 
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TABLE I: For each of the equilibrium polytropic stars used 
to test the code, the ratio of the ADM mass M to the TOV 
maximum mass M max xov, the ratio of the polar (i? po iar) to 
the equatorial (iiequat) coordinate radius, the ratio of the ro- 
tational kinetic energy T to the gravitational potential en- 
ergy \W\, and the ratio of the angular velocity on the equator 
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FIG. 3: The deviations from known equilibrium values of the 
density p and 4-metric ij) a g for various grids. As in Fig [5] 
Ldev(f) = L2(f — / ex )/L2(/ ex ). In the top panel, the metric 
is fixed and the fluid is evolved. In the other panels, the 
fluid is fixed and the metric is evolved. In the middle panel, 
the computational domain contains the stellar surface. In the 
bottom panel, the domain is entirely inside the star. 



C. Equilibrium Relativistic Stars 

In the previous tests, the matter flow did not affect the 
spacetime metric. For our next test, we evolve equilib- 
rium relativistic stars. Self-gravity is obviously of fun- 
damental importance for such objects, so the fluid and 
spacetime should now be evolved together, and we ex- 
pect significant feedback in both directions. For this 
test, we evolve isolated stars with polytropic equation 
of state P — Kop 1+1 / n , where we choose n = 1. We 
choose units such that G = c = Kq = 1. In these units, 
the maximum mass for a nonrotating n — 1 polytrope is 
AfmaxTOV = 0.164. This is the mass of a polytrope that 
has the critical central density p cr it = 0.32. Equilibrium 
configurations for rotating stars were computed using the 
code of Cook, Shapiro, and Teukolsky |67j |. 

We study four stars, whose properties are summarized 
in Table I. Stars A and B are nonrotating. Star A has 
central density (2/3)p cr it, and star B has central density 
(4/3)/? cr it- From the turning point theorem [68| . we infer 
that star A is stable and star B is unstable. These stars 



have octant symmetry and are evolved using fluid grids 
that cover one octant. Stars C and D are rotating, and we 
evolve them using fluid grids that cover the upper half- 
plane. Star C has the same rest mass as star A, but it 
rotates rigidly at an angular speed of 70% its mass shed- 
ding limit. From the turning point theorem, we know 
this star is secularly stable, and we expect it to be dy- 
namically stable as well. Star D is a differentially rotat- 
ing, hypermassive star similar to one studied by other 
groups [69l Frol ] . It is known from the simulations of these 
groups to be dynamically stable. 

First, we investigate the convergence of the fluid and 
PS codes, and in particular the effect of the discontinu- 
ity in the gradient of p on the convergence of the PS 
code. For this, it is sufficient to look at the spherically 
symmetric star A configuration. This star has a nonzero 
density at any point whose coordinate distance r from 
the star's center is less than the stellar radius R in these 
coordinates. 

We carry out three sets of convergence tests, with re- 
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FIG. 4: A neutron star with a realistic equation of state with 
both matter and metric evolved. Two runs are shown, corre- 
sponding to lower resolution PS and fluid grids and higher res- 
olution grids. We plot the normalized L2 norm of the change 
in density and the normalized L2 norm of the violation in the 
generalized harmonic constraints. 
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FIG. 5: The central density (relative to its initial value) for 
simulations of the four isolated stars. At the initial time, each 
star is subjected to a perturbation in the form of a one percent 
pressure depletion throughout the star. The inset shows the 
same thing but with a larger density range, so that later stages 
of the collapse of star B can be seen. 



suits shown in Fig [3] First, we hold the metric fixed and 
evolve the fluid using three resolutions of the fluid grid. 
As expected, we see second-order convergence. Because 
we plot the L2 norm of the deviation of p from its ini- 
tial value, the frequency of oscillations seen in the plot is 
twice the fundamental radial mode frequency. 

Next, we hold the hydrodynamic variables fixed and 
evolve the metric using the PS code. We fix the general- 
ized harmonic gauge source functions H a at their initial 
values, which are chosen to make the spacetime station- 
ary in our coordinates. To study the effects of the surface, 
we do two sets of runs. First, we evolve the metric on 
a domain that contains the stellar surface in the domain 
interior. The domain chosen is a spherical shell with in- 
ner radius at 0.7R and outer radius at 1.3i?. The middle 
panel of Fig [3] shows the results for PS grids with 8, 12, 
and 16 radial collocation points. The nonsmoothness of 
the metric at the surface reduces the convergence to sec- 
ond order. 

Finally, we evolve the metric in a PS domain in the 
interior of the star. The domain chosen is a cube at 
the center of the star with sides of length 0.86i?. We 
evolve with 8, 11, and 14 collocation points in each direc- 
tion. The resulting error is plotted in the bottom panel 
of Fig [3] The error should decrease exponentially as the 
number of points is increased, and we do see extremely 
rapid convergence. 

To test the full code, we next evolve a star in time, 
allowing both the hydrodynamic fields and the metric to 
change and influence each other. First, we check that 
the code still converges. To demonstrate the generality 
of the hydrodynamics code, we pick a more complicated 
equation of state for the star used in this test. (Of course, 
we have also checked that this test works with the poly- 
trope star A.) We choose the SLy neutron star equation 
of state 7l| to describe -P C oid, using the fitting function 



approximation presented in [23| . and we choose r t h = 2. 
The star has a central density of 10 15 g cm -3 and a mass 
of 1AM q. In Fig. 01 we show results from evolving this 
star at two resolutions. In the low resolution run, the 
fluid grid has 16 3 gridpoints covering one octant, while 
the PS grid is an ll 3 cube surrounded by 7 spherical 
shells, each with 6 radial collocation points and a maxi- 
mum angular I of 8. In the high resolution run, the fluid 
grid has 24 3 points, while for the PS grid we add one 
basis function in each direction in each domain. Fig. [4] 
shows the deviation of the density from its initial profile 
and the violation in the generalized harmonic constraint 
equations, the latter measured by the normalized con- 
straint energy ||C|| defined by Eq. (71) in [1^|). We see 
that both measures of error decrease significantly as res- 
olution is increased. (In this case, we have no expected 
convergence rate for comparison. However, we note that 
in Fig. [31 Ldcv(p) decreases as if it were converging to 
second-order in fluid grid spacing.) 

As a final test, we evolve all of the polytropic stars (A 
B C and D) for many dynamical times to check that our 
code can distinguish stable from unstable stars and can 
evolve stable stars accurately for long times. We evolve 
stars A and B using octant symmetry to restrict the fluid 
evolution to a single octant. For these stars, we divide 
the computational domain into 8 fluid domains with 18 3 
gridpoints each and 8 PS domains with roughly 10 3 collo- 
cation points on an average domain. For stars C and D, 
we utilize equatorial symmetry to restrict the fluid evolu- 
tion to the upper half-plane. We use 32 fluid domains of 
19 3 gridpoints and 32 PS domains of roughly 10 3 points 
each. For all four stars, the PS grid consists of a cube at 
the star's center, a spherical shell well outside the star at 
the outer boundary, and the region in between covered by 
several layers of cubed spheres, one of which contains the 
stellar surface. The generalized harmonic gauge source 
functions H a are again chosen to be fixed in time at their 
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initial values. In order to test the stability of each star, 
we perturb the initial data by reducing the pressure ev- 
erywhere by one percent. The induced violation of the 
constraint equations is small, so we do not adjust the 
metric to re-solve the constraints. 

In Figure 03 we plot the central density as a function 
of time for each star. As anticipated, stars A, C, and D 
are stable, while star B collapses. (We do not attempt 
to follow the collapse of star B to late times because 
we expect our choice of H a to be poor once the density 
profile changes drastically.) We see that the stable stars 
remain close to their initial states for many dynamical 
times. (For example, we have evolved the hypermassive 
star D for twenty central rotation periods.) Of the stable 
stars, star A has density oscillations with the highest 
amplitude, probably because it is very close to the TOV 
turning point. The oscillations have a period of 16.7, 
which is close to the expected period for linear radial 
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oscillations of ~ 7p c = 15.7. 
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FIG. 6: A representation of the choices of grid resolutions 
used in our reported runs. Nf is the effective total num- 
ber of grid points on the fluid finite difference grid, which in- 
cludes "mirror" points given by the reflection symmetry about 
the equatorial plane; Nps is the total number of collocation 
points on the spectral grid. 



IV. BLACK HOLE NEUTRON STAR BINARY 
INSPIRAL 

We now apply our code to model a black hole-neutron 
star (BHNS) binary, starting from the late stage of its 
inspiral. To date, there have been only a few success- 
ful simulations of these important systems in full general 
relativity [H, [H, HH, [36| . For this paper, we have chosen 
to evolve an equal mass system. Of course, we do not ex- 
pect such a system to occur frequently in nature, but we 
found this choice useful for two reasons. First, it allows 
us to make maximum use of our experience with binary 
black holes, which has largely focused on the equal mass 
case. Second, since some of the most important work 
on BHNS binaries with relativistic gravity (e.g. [3l|) has 
been done in the limit of extreme mass ratios, an equal 
mass system allows us to emphasize that we do not have 
this restriction. 

The code used to construct our initial data is described 
in (72| . For simplicity, and to facilitate comparisons with 
the results of other codes, we have modeled the neutron 
star as a r = 2 polytrope. At the start of the evolu- 
tion, the neutron star is chosen to be irrotational, and 
the black hole is chosen to have a zero quasi-local spin. 
(See [zl, [73| for details.) We use units such that the 
mass of the black hole and the neutron star are both 
one. (We define the mass of the star to be the ADM mass 
of an isolated star with the same baryonic mass as the 
one in the binary. Because the black hole has vanishing 
spin, we simply use the irreducible mass of the apparent 
horizon as the black hole mass.) At the initial time, we 
place the black hole and neutron star in nearly circular 
orbit with a binary separation of 24 in our coordinates 
and units. In particular, we place the center of rota- 
tion at the point (0,0,0), the maximum density inside 
the star at (—12,0,0), and the center of the horizon at 
(12.011, 0.097, 0). The rotation axis is the z-axis, and the 



initial period is P — 580. At this separation, the binary 
is expected to evolve a few orbits before the neutron star 
is disrupted. As we evolve, we use our dual frame sys- 
tem to hold the center of mass of the star (which nearly 
coincides with its maximum density point) at (—12, 0, 0) 
in moving coordinates, and we drive the horizon center 
to (12.011,0,0). After a quick adjustment, in which the 
black hole center moves onto y = 0, we find that the 
centers remain locked throughout the evolution. The co- 
ordinate mapping used to track the binary objects is a 
combination of two maps of type TZS described in Sec- 
tion [ITT3 One TZS, labeled TZS 1 , is chosen to have a 
rotation axis that runs through the origin, and its a and 
$ parameters are used to fix the x and y moving coor- 
dinates of the center of the apparent horizon. If the two 
binary objects were identical (as in the equal mass, non- 
spinning binary black hole problem), this one mapping 
would fix both objects by symmetry. In our case, how- 
ever, the black hole and neutron star are not identical, 
and asymmetries may develop. So we introduce a second 
mapping. The second TZS, labeled TZS 2 is chosen to have 
a rotation axis at x — 12.011, y — 0, the location of the 
black hole. This map does not affect the location of the 
black hole, but its a and $ are adjusted to fix the neutron 
star center of mass, defined as J x l Dd 3 x/ J Dd 3 x. In the 
case reported here, the second map is nearly an identity 
and so is not really needed. 

The fluid grid is chosen to be a cube of side length 16. 
(The diameter of the star is 11.5.) The pseudospectral 
grid consists of 7 cubes, 22 shells, 28 cylinders, and one 
cubed sphere. It extends to an outer radius of 326. In 
order to test the accuracy and convergence of the inspiral 
simulation, we vary the resolutions of both grids as shown 
in Fig. [51 We evolve with three fluid grids, having 48 3 x i, 
72 3 x |, 96 3 x \ total gridpoints. (The \ factors come 
from our use of equatorial symmetry to evolve only above 
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FIG. 7: The normalized constraint energy ||C|| defined by 
Eq. (71) in [39J. In this plot and all plots below, we use units 
in which the initial irreducible mass of the black hole is unity. 
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FIG. 8: The baryonic density and polytropic constant of the 
neutron star, normalized to their initial values. 



the equator.) We also use three PS grids, with 43 3 , 48 3 , 
and 54 3 collocation points. The 43 3 and the 48 3 grids 
(and also the 48 3 and 54 3 grids) differ by the addition 
of one collocation point (i.e. one new basis function) in 
each direction on each domain, with the exception that 
the azimuthal direction on cylinders is incremented by 
two points (which corresponds to adding a new sine and 
cosine basis function), as is needed to see convergence. 
To test boundary effects, the middle resolution (22 in 
Fig. [5]) is also run on a pseudospectral grid with outer 
radius of 400. 

Each resolution is run on 32 processors on the Caltech 
SHC cluster. For resolution 22, it takes 37 hours (1184 
CPU hours) to evolve one orbit at the initial separation. 
At resolution 33, it takes 100 hours (3200 CPU hours) to 
go one orbit. 

We evolve to t = 1000 (a little under two periods) at 
each resolution. The results are shown in figures El [51 
El and [TO] In Fig. [71 we plot the normalized constraint 
energy ||C|| as a function of time. We see that the con- 
straints converge quickly with PS resolution. Changing 
the fluid grid within the range studied has a very small 
effect on this diagnostic, i.e. the results for Runs 21 and 
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FIG. 9: The irreducible mass of the black hole. 
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FIG. 10: The proper separation between the black hole hori- 
zon and the neutron star center. 



23 (not plotted) would nearly coincide with that of 22. 
In the top panel of Fig. [8j we plot the central baryonic 
density normalized to its initial value. We find that the 
density converges with fluid grid resolution, but is insen- 
sitive to PS resolution over the range studied, i.e. Runs 
12 and 32 look like 22. At the highest fluid resolution, 
the central density is nearly constant, with only a small 
downward drift, over the roughly two orbits of evolution 
shown. Most of the unphysical density change is caused 
by spurious heating. This heating can be measured by 
the effective polytropic constant k = P/ p T . Since there 
are no physical shocks during this part of the evolution, k 
should be exactly constant. We plot k at the center of the 
star, normalized to its initial value, in the bottom panel 
of Fig[8j We find that the amount of spurious heating is 
very small, giving us confidence in our fluid evolutions. 

In FigEl we plot M- m: the irreducible mass of the black 
hole. Perhaps surprisingly, we find that the error in M; rr 
is sensitive to both grids, so to see convergence it is nec- 
essary to increase the resolutions of both grids simulta- 
neously. Therefore, we plot M- m for Runs 11, 22, and 33. 
We see that most of the increase in black hole mass con- 
verges away. In Fig 1101 we plot the proper separation d 
between the neutron star's center of mass and the point 
on the horizon on the a;-axis facing the star. Once again, 
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FIG. 11: The mass and angular momentum of the binary. 
Solid lines are ADM surface integrals. Dashed lines show the 
expected changes due to losses by the observed gravitational 
waves. 

errors from both grids contribute, so one must use Runs 
11, 22, and 33 to see convergence. We see that even our 
runs at higher resolution contain some eccentricity, pre- 
sumably an artifact of the approximations, particularly 
the choice of zero radial velocity, used to construct the 
initial data. (See [72| , in which we remove most of this 
eccentricity by adding an initial radial infall.) 

In Fig. [TTJ we plot the total ADM mass Madm and 
angular momentum Jadm of the system as a function of 
time, as measured by surface integrals carried out on the 
surface r — 300. We also monitor the gravitational wave 
emission via the Newman- Penrose scalar 1P4. (See [42| . 
for details on how ^4 is extracted from the evolution 
data.) After an initial burst of "junk" radiation, the true 
gravitational wave signal can be observed. Nearly all of 
the flux is carried in the I = 2, m = ±2 modes. The wave 
signal is a sinusoid whose amplitude and frequency are 
nearly constant over the course of our evolution. Reading 
off the amplitude, frequency, and time of arrival from the 
?/>4 measurements, we compute the predicted rate of mass 
and angular momentum loss. These are compared with 
the actual rates in Fig. [TTJ and we see that they agree 
well. The energy flux from a point mass binary with 
M = 2 and d = 24 is M = 1.61 x 10" 6 . The numerical 
flux is M = 1.55 x 10 -6 , so our numerical gravitational 
waves are reasonable. To test the effects of the outer 
boundary, we reran the 22 case with outer radius at 420 
and found no significant change in Madm, Jadm, or rtp4. 
Also, we have varied the extraction radii and found very 
small changes (apart from retardation effects) in Madm, 
Jadm, and rip^ for r between 150 and 400. 

V. BLACK HOLE NEUTRON STAR BINARY 
MERGER 

By continuing the evolution of the BHNS system, our 
code can simulate the merger of the two objects and 
determine the final stationary state. However, we find 



that we must introduce several modifications to our evo- 
lution techniques during the merger phase, which we 
find begins after about four orbits of inspiral, at around 

t = tdisr = 1700. 

The most obvious adjustment to be made is in the ex- 
tents of the fluid grid. When the neutron star starts to 
disrupt, it is no longer possible to confine this grid to a 
small box around the star. So we must regrid at certain 
times to keep the matter on the grid. Actually, regrid- 
ding may also be necessary from time to time during long 
inspirals, because the neutron star grows in the moving 
coordinates as the scaling parameter a in the 1ZS 1 map- 
ping decreases. In our simulation, we need to expand 
the fluid grid three times, at t = 1080, t — 1540, and 
t = tdisr- The final fluid grid must cover a region around 
the neutron star and the black hole. Thus, some of the 
savings of our two-grid approach is lost. Even so, the 
fluid grid does not need to extend all the way into the 
wave zone, as it does in pure finite difference codes. 

Crucial adjustments must also be made to the dual 
frame control system. During the disruption, it is no 
longer appropriate to fix the center of the neutron star. 
In a sense, the BHNS problem is easier for our code than 
the binary black hole problem. In the binary black hole 
case, the presence of two excision zones forces us to fix 
the locations of both holes in moving coordinates, leading 
to an very distorted moving coordinate system as the two 
holes approach each other in inertial coordinates. In the 
BHNS case, we can "release" the neutron star as it is 
disrupted, and allow it to fall into the black hole. This 
is done by setting $(t) = <I>(tdisr) and a(t) = a(tdisr) for 
t > tdisr in the 1ZS 2 mapping used to fix the star. 

Because we use excision, we must continue to fix the 
black hole's location on our grid. However, if we were 
to try do do this with the 1ZS 1 mapping, the parameter 
a would approach zero as the hole moved towards the 
origin. Instead, we switch our coordinate control system 
at i — tdisr- We cease evolving the 1ZS 1 parameters, 
fixing a(t) = a(t disr ) and $(t) = $(t dis r) + $ (tdisr )(t - 
tdisr) for t > tdisr • (Letting the time derivative of <f> 
be continuous makes the coordinate evolution smoother.) 
We then compose this mapping with another, a simple 
translation: x l — x % + C % . C x and C v are used to fix the 
center of the horizon. 

In addition to moving, the horizon changes in size and 
shape, becoming much larger as it accretes matter and 
becoming very distorted during the merger. We find that 
our code becomes less accurate and less stable the deeper 
the PS grid extends inside the hole. Therefore, it is im- 
portant to place the excision boundary as close to the 
apparent horizon as possible. We do this by introducing 
another coordinate mapping that controls the size and 
shape of the horizon. This mapping has the form 

r = r-f(r)Y l \ lm Y lm (6,<p) , (20) 

lm 

where r is the coordinate distance from the horizon center 
(12.011,0,0), Xim are functions of time, and f(r) is a 




FIG. 12: Snapshots of the density and horizon shape on the 
equatorial plane at times t = 1730, t = 1850, and t = 1970. 
Shades of grey represent p/p ma x, the density relative to its 
current maximum. The black object is the apparent horizon. 
The density and horizon are shown in inertial coordinates. 
The distorted rectangle is the outer edge of the fluid grid, 
which is a fixed rectangle in moving coordinates but moves in 
the inertial frame. 



smooth function of r. (We use a Gaussian.) The A; m 
parameters are used to drive the corresponding moment 
of the horizon's shape in moving coordinates. For the 
simulation described below, the sum in Eq. I|20p extends 
from I = to I = 6 modes, omitting I — 1 modes, which 
are controlled by the translation. At t = tdisr, we set 
A; TO = for all I > 0. A o is set to a value that doubles the 
diameter of the horizon in moving coordinates, effectively 
increasing the grid coverage near the black hole. Because 
this is a discontinuous change of coordinates, we must 
interpolate the evolution variables onto the new grid at 
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FIG. 13: The evolution of the neutron star and black hole 
masses from the beginning of the star's disruption to the final 
stationary state. Here, Mo is the rest mass of the neutron star 
matter outside the hole, M\ TI is the black hole's irreducible 
mass, and M c is the Christodoulou mass of the black hole, 
defined as *J M? r + J 2 /(4M^ r ), where J is the black hole's 
spin. 
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FIG. 14: The trajectory in inertial coordinates of the neutron 
star center from t = to t = tdisn and that of the black 
hole horizon center from t = to t = 3000. The black hole 
line switches from dashed to solid at t — £disr- Both the 
black hole and the neutron star remain on the equatorial plane 
throughout the evolution. 



^ — ^disr- 

Finally, we find minor improvements by adjusting the 
generalized harmonic gauge functions H a during the 
merger. In the future, we hope to do this using dynamical 
gauge conditions. For now, we simply damp them to zero 
by setting H a (t) = H a (Q) e - ( ' < - t " tdi ^/ T ^ , where r = 100. 
Our evolutions seem to be fairly insensitive to r; it can 
be varied by an order of magnitude with no significant 
effects. 

For our merger simulation, we begin at t = 1000 with 
the inspiral data computed on the grids corresponding 
to resolution 22 and PS outer radius of 400. We then 
continue the evolution of this system to later times. The 
evolution from t — 1000 to t — 1700 is very similar to the 



13 




1000 



2000 



3000 



FIG. 15: The gravitational wave signal for the entire inspiral 
and merger (excluding the initial burst of "junk" radiation). 
The quadrupole contribution dominates, so we plot only the 
I — m — 2 amplitudes. 



evolution from t — to t = 1000 described in Section HVl 
At t = 1700, matter is starting to flow off the neutron 
star, so we introduce the changes described above, con- 
tinuing the evolution on 80 processors in order to ac- 
commodate the larger fluid grid needed. The resolution 
on the final fluid grid is Ax = 0.4, so that there are 20 
fluid grid points across the diameter of the excision zone. 
On the LONI Queen Bee cluster, the code must run for 
43 hours (3440 CPU hours) to evolve from t = 1700 to 
t = 2000, by which time most of the matter has fallen 
into the hole. 

Three snapshots of the neutron star density are shown 
in Fig. [T3J A dense stream of matter flows from the star 
to the black hole, so that there quickly come to be two 
peaks in the matter density: one at the center of the 
neutron star and the other at the point where the stream 
reaches the hole. The neutron star continues to fall closer 
to the hole, even as it rapidly loses mass through the 
matter stream. The density peak corresponding to the 
neutron star core disappears at around t — 1800. At 
around t = 1850, the matter stream starts to close in a 
ring around the black hole. A shock forms at this time 
where matter flowing towards the hole intersects matter 
flowing around the hole. This can be seen in the sharp in- 
ner edge of the grey swath in the middle panel of Fig. [T2J 
Matter falls rapidly into the hole from t = 1750 until 
t = 2000. 

In Fig. [13l we plot the baryon rest mass Mq on the grid 
together with the apparent horizon irreducible mass M; rr 
as functions of time starting from t = tdisr- We find that 
virtually all of the matter falls promptly into the black 
hole; Mq drops from its initial value near unity until it 
stabilizes about 2 x 10~ 4 , indicating a very small post- 
merger accretion disk. This result is consistent with the 
small disks found by other groups (particularly Etienne et 
al [36[ , who have also simulated the equal mass case, but 
also the most recent results of Yamamoto, Shibata, and 
Taniguchi [35]). At t = 2200, the matter is so sparse 
that we stop evolving the fluid and continue the metric 
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FIG. 16: The constraint energy ||C||, black hole mass M c , 
and baryonic rest mass Mo for the merger phase in evolutions 
carried out at two different resolutions. 



evolution using our vacuum Einstein code. 

The final M- m of the black hole is very close to 2, and 
its final spin is J/M 2 ss 0.7. (The computations used to 
obtain the black hole spin are described in the appen- 
dicies of [z3|.) From M m and J, we can compute the 
Christodoulou mass, and it is found to be M c « 2.15. 
That M c is greater than 2 is presumably a result of the 
error in our calculation and is consistent with the ob- 
served constraint violation (see below). The complete 
trajectories in inertial coordinates of the black hole hori- 
zon center and the neutron star center of mass are plotted 
in Fig. [T3J In Fig. Qj)J we plot the I = m = 2 component 
of rip4 for the complete evolution, extracting the wave at 
r = 350. 

To test the stability and convergence rate of our merger 
algorithm, we have run the merger phase 1700 < t < 2050 
at a higher resolution. For this run, we increase the num- 
ber of points on the fluid grid by 40% and the number of 
PS collocation points by one per axis per domain, keep- 
ing the extents of both grids fixed. The code runs about 
four times slower on the larger grid. The effects of this 
grid change are shown in Fig. 1161 The constraint vio- 
lation peaks at t ss 1900 at about ||C|| =0.1 for the 
low resolution run and at around ||C|| = 0.07 for the 
high resolution run. The constraint violations decline 
thereafter. We note that, since the constraint energy is 
strongly peaked near the black hole, these numbers de- 
pend sensitively on how the norm is taken. The quantity 

1/2 

||C|| is an L2-type norm of the form (/ f 2 dV) , where 
/ represents the generalized harmonic constraints. (See 
Eq. 53 and Eq. 71 of 39].) If we instead use an Ll- 
type norm of the form J \ f\ dV , we find that this norm, 
when appropriately normalized, peaks at 0.02 for the low 
resolution run and 0.01 for the high resolution run. The 
Hamiltonian and momentum constraints (integrated with 
Ll-type norms) show similar levels of violation. The rel- 
ative momentum constraint violation peaks at about 0.02 
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for the low resolution run and 0.01 for the high resolu- 
tion run. The relative Hamiltonian constraint violation is 
about 0.02 for both resolutions — convervence is not seen 
because this error is present in the data at t = i^isr from 
which both runs start. Of course, these numbers are also 
very sensitive to the way the constraints are summed and 
normalized. The post-merger black hole mass is slightly 
lower at higher resolution. The accretion rate is also 
lower, but only after about 98% of the rest mass has al- 
ready fallen into the hole. In neither run do we find a 
long-lived massive disk. 

VI. CONCLUSIONS 

We have shown that our code can evolve inspiraling 
BHNS binaries with high accuracy at fairly low compu- 
tational cost. This enables us to begin simulations at 
relatively large binary separations. Long inspirals may 
turn out to be very important for accurately modeling 
mergers; merger simulations by the Illinois_group show 
some sensitivity to the initial separation [36], with the 
implication that starting too close to the merger can lead 
to underestimating the mass of the post-merger accre- 
tion disk. We have also shown that we can simulate the 
merger of these binaries, although our accuracy of these 
simulations is not as good so far than that of our inspirals. 
Also, there is every reason to believe that the techniques 
described here would work just as well for binary neutron 
stars. 

An important next step is to demonstrate that our 
code can evolve more general BHNS binaries. We have 



recently begun evolving such systems with different mass 
ratios and black hole spins, and so far we have had little 
difficulty in evolving the inspirals. We hope to report 
on these simulations in the near future. Also, we plan 
to study ways to improve the accuracy of our merger 
simulations. 

It has only been two years since the first fully relativis- 
tic BHNS merger simulations were reported, and as yet 
very little of the parameter space has been studied. In 
particular, there is a pressing need to study the effects 
of the black hole spin and the neutron star equation of 
state. There are preliminary indications (e.g. [§, [34|) 
that both of these have important effects on the post- 
merger accretion disk mass and the gravitational wave 
signal. Also, both of the other groups currently perform- 
ing BHNS merger simulations use very similar techniques 
(BSSN, moving punctures). It will be useful to compare 
their results with those obtained using very different tech- 
niques, like the ones reported here. 
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